// Calculations of Berry curvature in two different ways

#include "tbparameters.hpp"
cx_mat berry_sumrule(cx_mat Vx, cx_mat Vy, imat mask, vec E){
  
  int fdim = E.n_elem;
  cx_mat Rx(fdim,fdim);
  cx_mat Ry(fdim,fdim);

  cx_double cmplxi = cx_double(0.0,1.0);

  Rx.zeros();
  Ry.zeros();

  for (int i = 0; i < fdim; i++)
    for (int j = 0; j < i; j++){
      //cout << i << " " << j << " " << mask.n_rows << " " << fdim << "\n";
      if (mask(j,i) == 0 ){
	Rx(j,i) = -cmplxi*Vx(j,i)/(E(j)-E(i));
	Ry(j,i) = -cmplxi*Vy(j,i)/(E(j)-E(i));
	Rx(i,j) = conj(Rx(j,i));
	Ry(i,j) = conj(Ry(j,i));
    }
    }

  cx_mat B(fdim,fdim);
  // Rx Ry - Vy Vx
  B = cmplxi*(Rx * Ry - Ry * Rx);
  B = (B+B.t())/2;
  
  B = B % mask;

  return B;


}

cx_mat berry(cx_mat *eigvec, imat mask, int nplaq,double plaqarea, cx_mat L[]){
  
  //cx_mat L[nplaq];
  vec eigval[nplaq];
  int i,j;

  int fdim=eigvec[0].n_rows;
  int rdim=fdim;
  cx_mat v[nplaq];

  mat t;

 
  cx_mat y1,y2;
  vec s;

  cx_mat eigvec_dag[nplaq];

  for (i=0;i<nplaq;i++)
    eigvec_dag[i]=eigvec[i].t();


  for (i=0;i<nplaq;i++){
    if (i < nplaq-1 )
      j=i+1;
    else
      j=0;

    L[i] = eigvec_dag[i] * eigvec[j];
    L[i] = L[i] % mask;
    svd(y1,s,y2,L[i]);
    L[i] = y1 * y2.t();
  }

  //cout << abs(L[0].submat(4,4,7,7))<<"\n";

  cx_mat B=L[0] % mask;
  for (i=1;i<nplaq;i++)
    B = B*(L[i] % mask);


  // 1-iB = U
  cx_mat Id;

  Id.eye(fdim,fdim);
  B = cx_double(0.0,1.0)*(Id - B);


  // Below the division by (pi*pi/4) is due to the fact that 
  // k being used on the grid is related to the argument of 
  // plane wavees for tight binding by (k*pi/2). See Vogl's paper
  B = (B+B.t())/2 / plaqarea / (pi*pi/4);
  


  //cout << abs(B.submat(4,4,7,7))<<" B\n";

  return B;
}

